# PoP - Policing Socio-Geographic Boundaries and Inequality
# script for creating Figure 1: Influence of Racial Boundary on Logged Arrests (Standardized). 

# load packages
suppressPackageStartupMessages(
  
  {
    library(dplyr)
    library(tidyverse)
    library(ggplot2)
    library(sensemakr)
    library(haven)
    library(readxl)
    library(readr)
    library(areal)
    library(car)
    library(estimatr)
    library(magrittr)
    library(texreg)
    library(sandwich)
    library(jtools)
    library(ggthemes)
    library(meta)
  }
)

# read in data by city 
# Load Atlanta final dataset
load("atl_final.RData")

# prepare final variables by city

# log and +1 to variables (pop already logged)
atl_fin = atl_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
atl_fin$larrest_sd <- atl_fin$larrests - (mean(atl_fin$larrests)/sd(atl_fin$larrests))
atl_fin$lmisdemeanors_sd <- atl_fin$lmisdemeanors - (mean(atl_fin$lmisdemeanors)/sd(atl_fin$lmisdemeanors))
atl_fin$lfelonies_sd <- atl_fin$lfelonies - (mean(atl_fin$lfelonies)/sd(atl_fin$lfelonies))
atl_fin$lnonviolent_sd <- atl_fin$lnonviolent - (mean(atl_fin$lnonviolent)/sd(atl_fin$lnonviolent))
atl_fin$lviolent_sd <- atl_fin$lviolent - (mean(atl_fin$lviolent)/sd(atl_fin$lviolent))
atl_fin$lsociety_sd <- atl_fin$lsociety - (mean(atl_fin$lsociety)/sd(atl_fin$lsociety))
atl_fin$lperson_sd <- atl_fin$lperson - (mean(atl_fin$lperson)/sd(atl_fin$lperson))
atl_fin$lproperty_sd <- atl_fin$lproperty - (mean(atl_fin$lproperty)/sd(atl_fin$lproperty))

# rename racial and economic boundary variable
atl_fin$atl_white_blv <- atl_fin$p_race_white_blv
atl_fin$atl_ses_blv <- atl_fin$ses_blv


# Load Austin final dataset
load("aus_final.RData")


# log and +1 to variables (pop already logged)
aus_fin = aus_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
aus_fin$larrest_sd <- aus_fin$larrests - (mean(aus_fin$larrests)/sd(aus_fin$larrests))
aus_fin$lmisdemeanors_sd <- aus_fin$lmisdemeanors - (mean(aus_fin$lmisdemeanors)/sd(aus_fin$lmisdemeanors))
aus_fin$lfelonies_sd <- aus_fin$lfelonies - (mean(aus_fin$lfelonies)/sd(aus_fin$lfelonies))
aus_fin$lnonviolent_sd <- aus_fin$lnonviolent - (mean(aus_fin$lnonviolent)/sd(aus_fin$lnonviolent))
aus_fin$lviolent_sd <- aus_fin$lviolent - (mean(aus_fin$lviolent)/sd(aus_fin$lviolent))
aus_fin$lsociety_sd <- aus_fin$lsociety - (mean(aus_fin$lsociety)/sd(aus_fin$lsociety))
aus_fin$lperson_sd <- aus_fin$lperson - (mean(aus_fin$lperson)/sd(aus_fin$lperson))
aus_fin$lproperty_sd <- aus_fin$lproperty - (mean(aus_fin$lproperty)/sd(aus_fin$lproperty))

# rename racial and economic boundary variable
aus_fin$aus_white_blv <- aus_fin$p_race_white_blv
aus_fin$aus_ses_blv <- aus_fin$ses_blv


## Load Boston data
load("bos_final.RData")


# log and +1 to variables (pop already logged)
bos_fin = bos_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime + 1))



# create scaled DVs (to account for different pop. sizes and density in blocks)
bos_fin$larrests_sd <- bos_fin$larrests - (mean(bos_fin$larrests)/sd(bos_fin$larrests))
bos_fin$lmisdemeanors_sd <- bos_fin$lmisdemeanors - (mean(bos_fin$lmisdemeanors)/sd(bos_fin$lmisdemeanors))
bos_fin$lfelonies_sd <- bos_fin$lfelonies - (mean(bos_fin$lfelonies)/sd(bos_fin$lfelonies))
bos_fin$lnonviolent_sd <- bos_fin$lnonviolent - (mean(bos_fin$lnonviolent)/sd(bos_fin$lnonviolent))
bos_fin$lviolent_sd <- bos_fin$lviolent - (mean(bos_fin$lviolent)/sd(bos_fin$lviolent))
bos_fin$lsociety_sd <- bos_fin$lsociety - (mean(bos_fin$lsociety)/sd(bos_fin$lsociety))
bos_fin$lperson_sd <- bos_fin$lperson - (mean(bos_fin$lperson)/sd(bos_fin$lperson))
bos_fin$lproperty_sd <- bos_fin$lproperty - (mean(bos_fin$lproperty)/sd(bos_fin$lproperty))

# rename racial and economic boundary variable
bos_fin$bos_white_blv <- bos_fin$p_race_white_blv
bos_fin$bos_ses_blv <- bos_fin$ses_blv

## Load Chicago data
load("chi_final.RData")

# log and +1 to variables (pop already logged)
chi_fin = chi_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(property_crime+1),
         lviolentcrime = log(violent_crime+1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
chi_fin$larrest_sd <- chi_fin$larrests - (mean(chi_fin$larrests)/sd(chi_fin$larrests))
chi_fin$lmisdemeanors_sd <- chi_fin$lmisdemeanors - (mean(chi_fin$lmisdemeanors)/sd(chi_fin$lmisdemeanors))
chi_fin$lnonviolent_sd <- chi_fin$lnonviolent - (mean(chi_fin$lnonviolent)/sd(chi_fin$lnonviolent))
chi_fin$lsociety_sd <- chi_fin$lsociety - (mean(chi_fin$lsociety)/sd(chi_fin$lsociety))
chi_fin$lfelonies_sd <- chi_fin$lfelonies - (mean(chi_fin$lfelonies)/sd(chi_fin$lfelonies))
chi_fin$lperson_sd <- chi_fin$lperson - (mean(chi_fin$lperson)/sd(chi_fin$lperson))
chi_fin$lproperty_sd <- chi_fin$lproperty - (mean(chi_fin$lproperty)/sd(chi_fin$lproperty))
chi_fin$lviolent_sd <- chi_fin$lviolent - (mean(chi_fin$lviolent)/sd(chi_fin$lviolent))


# rename racial and economic boundary variable
chi_fin$chi_white_blv <- chi_fin$p_race_white_blv
chi_fin$chi_ses_blv <- chi_fin$ses_blv

## Load Louisville data
load("lou_final.RData")


# log and +1 to variables (pop already logged) (no misdemeanors vs felonies for louisville)
lou_fin = lou_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property+1),
         lviolentcrime = log(crime_violent+1),
         lviolent = log(violent_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
lou_fin$larrest_sd <- lou_fin$larrests - (mean(lou_fin$larrests)/sd(lou_fin$larrests))
lou_fin$lnonviolent_sd <- lou_fin$lnonviolent - (mean(lou_fin$lnonviolent)/sd(lou_fin$lnonviolent))
lou_fin$lsociety_sd <- lou_fin$lsociety - (mean(lou_fin$lsociety)/sd(lou_fin$lsociety))
lou_fin$lviolent_sd <- lou_fin$lviolent - (mean(lou_fin$lviolent)/sd(lou_fin$lviolent))
lou_fin$lperson_sd <- lou_fin$lperson - (mean(lou_fin$lperson)/sd(lou_fin$lperson))
lou_fin$lproperty_sd <- lou_fin$lproperty - (mean(lou_fin$lproperty)/sd(lou_fin$lproperty))

# rename racial and economic boundary variable
lou_fin$lou_white_blv <- lou_fin$p_race_white_blv
lou_fin$lou_ses_blv <- lou_fin$ses_blv


## Load Milwaukee data
load("mil_final.RData")


# log and +1 to variables (pop already logged)
mil_fin = mil_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
mil_fin$larrest_sd <- mil_fin$larrests - (mean(mil_fin$larrests)/sd(mil_fin$larrests))
mil_fin$lmisdemeanors_sd <- mil_fin$lmisdemeanors - (mean(mil_fin$lmisdemeanors)/sd(mil_fin$lmisdemeanors))
mil_fin$lfelonies_sd <- mil_fin$lfelonies - (mean(mil_fin$lfelonies)/sd(mil_fin$lfelonies))
mil_fin$lnonviolent_sd <- mil_fin$lnonviolent - (mean(mil_fin$lnonviolent)/sd(mil_fin$lnonviolent))
mil_fin$lviolent_sd <- mil_fin$lviolent - (mean(mil_fin$lviolent)/sd(mil_fin$lviolent))
mil_fin$lsociety_sd <- mil_fin$lsociety - (mean(mil_fin$lsociety)/sd(mil_fin$lsociety))
mil_fin$lperson_sd <- mil_fin$lperson - (mean(mil_fin$lperson)/sd(mil_fin$lperson))
mil_fin$lproperty_sd <- mil_fin$lproperty - (mean(mil_fin$lproperty)/sd(mil_fin$lproperty))

# rename racial and economic boundary variable
mil_fin$mil_white_blv <- mil_fin$p_race_white_blv
mil_fin$mil_ses_blv <- mil_fin$ses_blv


## Load Seattle data
load("sea_final.RData")

# log and +1 to variables (pop already logged)
sea_fin = sea_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
sea_fin$larrest_sd <- sea_fin$larrests - (mean(sea_fin$larrests)/sd(sea_fin$larrests))
sea_fin$lmisdemeanors_sd <- sea_fin$lmisdemeanors - (mean(sea_fin$lmisdemeanors)/sd(sea_fin$lmisdemeanors))
sea_fin$lfelonies_sd <- sea_fin$lfelonies - (mean(sea_fin$lfelonies)/sd(sea_fin$lfelonies))
sea_fin$lnonviolent_sd <- sea_fin$lnonviolent - (mean(sea_fin$lnonviolent)/sd(sea_fin$lnonviolent))
sea_fin$lviolent_sd <- sea_fin$lviolent - (mean(sea_fin$lviolent)/sd(sea_fin$lviolent))
sea_fin$lsociety_sd <- sea_fin$lsociety - (mean(sea_fin$lsociety)/sd(sea_fin$lsociety))
sea_fin$lperson_sd <- sea_fin$lperson - (mean(sea_fin$lperson)/sd(sea_fin$lperson))
sea_fin$lproperty_sd <- sea_fin$lproperty - (mean(sea_fin$lproperty)/sd(sea_fin$lproperty))

# rename racial and economic boundary variable
sea_fin$sea_white_blv <- sea_fin$p_race_white_blv
sea_fin$sea_ses_blv <- sea_fin$ses_blv


# Sensitivity Analysis with DV as Misdemeanor arrests
# atlanta
atl.model.2 = lm_robust(lmisdemeanors_sd ~ atl_white_blv + atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                          hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                          lviolentcrime + phys_bound + com_density, 
                        data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)


atl.model.2_sense = lm(lmisdemeanors_sd ~ atl_white_blv + atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                         hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                         lviolentcrime + phys_bound + com_density, 
                       data = atl_fin, subset = pop > 0)


# create df of atlanta coefficients
coef.atl.2 <- tidy(atl.model.2)

fit_cis_95 <- confint(atl.model.2, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(atl.model.2, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.atl.2 <- bind_cols(coef.atl.2, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef
atl_boundary_white <- coef.atl.2[1,]


# austin
aus.model.2 = lm_robust(lmisdemeanors_sd ~ aus_white_blv + aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                          hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                          lviolentcrime + phys_bound + com_density, 
                        data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)


aus.model.2_sense = lm(lmisdemeanors_sd ~ aus_white_blv + aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                         hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                         lviolentcrime + phys_bound + com_density, 
                       data = aus_fin, subset = pop > 0)


# create df of austin coefficients
coef.aus.2 <- tidy(aus.model.2)

fit_cis_95 <- confint(aus.model.2, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(aus.model.2, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.aus.2 <- bind_cols(coef.aus.2, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef
aus_boundary_white <- coef.aus.2[1,]


# boston
bos.model.2 = lm_robust(lmisdemeanors_sd ~ bos_white_blv + bos_ses_blv + lpop + p_race_white + age_15_35_male + 
                          hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + 
                          phys_bound + com_density, 
                        data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)


bos.model.2_sense = lm(lmisdemeanors_sd ~ bos_white_blv + bos_ses_blv + lpop + p_race_white + age_15_35_male + 
                         hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + 
                         phys_bound + com_density, 
                       data = bos_fin, subset = pop > 0)


# create df of bostin coefficients
coef.bos.2 <- tidy(bos.model.2)

fit_cis_95 <- confint(bos.model.2, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(bos.model.2, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.bos.2 <- bind_cols(coef.bos.2, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef
bos_boundary_white <- coef.bos.2[1,]

# chicago
chi.model.2 = lm_robust(lmisdemeanors_sd ~ chi_white_blv + chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                          hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                          lviolentcrime + phys_bound + com_density, 
                        data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi.model.2_sense = lm(lmisdemeanors_sd ~ chi_white_blv + chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                         hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                         lviolentcrime + phys_bound + com_density, 
                       data = chi_fin, subset = pop > 0)

# create df of chicago coefficients
coef.chi.2 <- tidy(chi.model.2)

fit_cis_95 <- confint(chi.model.2, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(chi.model.2, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.chi.2 <- bind_cols(coef.chi.2, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef
chi_boundary_white <- coef.chi.2[1,]

# louisville - does not have misdemeanor vs felony distinction

# milwaukee
mil.model.2 = lm_robust(lmisdemeanors_sd ~ mil_white_blv + mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                          hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                          lviolentcrime + phys_bound + com_density, 
                        data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

mil.model.2_sense = lm(lmisdemeanors_sd ~ mil_white_blv + mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                         hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                         lviolentcrime + phys_bound + com_density, 
                       data = mil_fin, subset = pop > 0)

# create df of milwaukee coefficients
coef.mil.2 <- tidy(mil.model.2)

fit_cis_95 <- confint(mil.model.2, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(mil.model.2, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.mil.2 <- bind_cols(coef.mil.2, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef
mil_boundary_white <- coef.mil.2[1,]


# seattle
sea.model.2 = lm_robust(lmisdemeanors_sd ~ sea_white_blv + sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                          hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                          lviolentcrime + phys_bound + com_density, 
                        data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

sea.model.2_sense = lm(lmisdemeanors_sd ~ sea_white_blv + sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                         hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                         lviolentcrime + phys_bound + com_density, 
                       data = sea_fin, subset = pop > 0)


# create df of seattle coefficients
coef.sea.2 <- tidy(sea.model.2)

fit_cis_95 <- confint(sea.model.2, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(sea.model.2, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.sea.2 <- bind_cols(coef.sea.2, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef
sea_boundary_white <- coef.sea.2[1,]


## combine 6 dfs (lousiville not in this analysis) w white boundary coefs
combined.coeffs <- rbind(atl_boundary_white, aus_boundary_white, bos_boundary_white, chi_boundary_white,
                         mil_boundary_white, sea_boundary_white)


combined.coeffs <- combined.coeffs %>% dplyr::select(-c(SE,statistic,p.value, conf.low, conf.high))

## make xaxis labels
Plot1labels <- c("Atlanta", "Austin", "Boston","Chicago","Milwaukee","Seattle")

## sensitivity analysis

# robustness value exercise 

sense_out1 = 
  sensemakr(atl.model.2_sense,
            treatment = "atl_white_blv",
            benchmark_covariates = "lviolentcrime",
            kd = 1)


sense_out2 = 
  sensemakr(aus.model.2_sense,
            treatment = "aus_white_blv",
            benchmark_covariates = "lviolentcrime",
            kd = 1)


sense_out3 = 
  sensemakr(bos.model.2_sense,
            treatment = "bos_white_blv",
            benchmark_covariates = "lcrime",
            kd = 1)


sense_out4 = 
  sensemakr(chi.model.2_sense,
            treatment = "chi_white_blv",
            benchmark_covariates = "lviolentcrime",
            kd = 1)

sense_out5 = 
  sensemakr(mil.model.2_sense,
            treatment = "mil_white_blv",
            benchmark_covariates = "lviolentcrime",
            kd = 1)


sense_out6 = 
  sensemakr(sea.model.2_sense,
            treatment = "sea_white_blv",
            benchmark_covariates = "lviolentcrime",
            kd = 1)

rvdf = data.frame(
  
  rv_val = c(as.numeric(sense_out1$sensitivity_stats$rv_q),
             as.numeric(sense_out2$sensitivity_stats$rv_q),
             as.numeric(sense_out3$sensitivity_stats$rv_q),
             as.numeric(sense_out4$sensitivity_stats$rv_q),
             as.numeric(sense_out5$sensitivity_stats$rv_q),
             as.numeric(sense_out6$sensitivity_stats$rv_q)),
  
  variable = combined.coeffs$Variable,
  
  bound = c("5x Violent \nCrime", "3x Violent \nCrime", "2x Crime", "19x Violent \nCrime", 
            "4220x Violent \nCrime","5x Violent \nCrime")
  
)


combined.coeffs$se = (combined.coeffs$conf.high_95 - combined.coeffs$Coefficient) / 1.96

out_meta = metagen(
  TE = combined.coeffs$Coefficient,
  seTE = combined.coeffs$se
)

# Atlanta: 5x violent crime
# Austin: 3x violent crime
# Boston: 2x crime
# Chicago: 19x violent crime
# Louisville: NA
# Milwaukee: 4220x violent crime
# Seattle: 5x violent crime


df_meta = data.frame(
  Variable = "Meta-analysis",
  Coefficient = out_meta$TE.random,
  df = out_meta$df.Q,
  outcome = "larrest_sd",
  se = out_meta$seTE.random,
  conf.low_95 = out_meta$TE.random - 1.96 * out_meta$seTE.random,
  conf.high_95 = out_meta$TE.random + 1.96 * out_meta$seTE.random,
  conf.low_90 = out_meta$TE.random - 1.645 * out_meta$seTE.random,
  conf.high_90 = out_meta$TE.random + 1.645 * out_meta$seTE.random
)

combined.coeffs = bind_rows(combined.coeffs, df_meta)
combined.coeffs$Variable = 
  factor(combined.coeffs$Variable, levels = c("Meta-analysis", "atl_white_blv","aus_white_blv","bos_white_blv",
                                              "chi_white_blv","mil_white_blv","sea_white_blv"))




misd_white_blv_plot_all_cities <- 
  ggplot(combined.coeffs, 
         aes(x = Variable, y = Coefficient)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 1) +
  geom_point(aes(x = Variable, 
                 y = Coefficient)) + 
  geom_errorbar(aes(x = Variable,
                    ymin = conf.low_90,
                    ymax = conf.high_90),
                width = 0,
                size = 1) +
  geom_errorbar(aes(x = Variable,
                    ymin = conf.low_95,
                    ymax = conf.high_95),
                width = 0,
                size = 0.6) +
  annotate("text",
           label =  paste0("RV: ", round(rvdf$rv_val[1:6], 2),
                           "\n", rvdf$bound[1:6]),
           y = combined.coeffs$Coefficient[1:6],
           x = c(2.6, 3.6, 4.5, 5.6, 6.7, 7.6),
           size = 3.75,
           family = "serif") +
  geom_hline(yintercept = combined.coeffs$Coefficient[7], 
             colour = gray(1/2), lty = 2,
             alpha = .4,
             size = .2) + ylim(-0.3, 0.8) + 
  ggtitle("Influence of White Racial Boundaries on Misdemeanor Arrests") + 
  xlab("Cities") + ylab("Racial Boundary Coefficients") + 
  scale_color_grey(start = .6, end = 0) + 
  theme_tufte(base_size = 14) + 
  scale_x_discrete(labels = c("Meta-Analysis",Plot1labels))

# make meta-analytic coefficient variable red
misd_white_blv_plot_all_cities <- misd_white_blv_plot_all_cities + 
  geom_point(data=subset(combined.coeffs, Variable == "Meta-analysis"), colour="red") +
  geom_errorbar(data = subset(combined.coeffs, Variable == "Meta-analysis"),
                mapping = aes(x = Variable,
                              ymin = conf.low_90,
                              ymax = conf.high_90),
                color="red",
                width = 0,
                size = 0.8) +
  geom_errorbar(data = subset(combined.coeffs, Variable == "Meta-analysis"),
                mapping = aes(x = Variable,
                              ymin = conf.low_95,
                              ymax = conf.high_95),
                color="red",
                width = 0,
                size = .4)

# expand size of plot so full labels can be seen
misd_white_blv_plot_all_cities <- misd_white_blv_plot_all_cities + 
  expand_limits(x= c(length(levels(combined.coeffs$Variable)) + 1))



ggsave(plot = misd_white_blv_plot_all_cities, filename = "misd_white_blv_plot_sens.png", 
       width = 8, height = 4, bg = 'white')



# Felony arrests
# atlanta
atl.model.3 = lm_robust(lfelonies_sd ~ atl_white_blv + atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                          hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                          lviolentcrime + phys_bound + com_density, 
                        data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)


atl.model.3_sense = lm(lfelonies_sd ~ atl_white_blv + atl_ses_blv + lpop + p_race_white + age_15_35_male + 
                         hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                         lviolentcrime + phys_bound + com_density, 
                       data = atl_fin, subset = pop > 0)


# create df of atlanta coefficients
coef.atl.3 <- tidy(atl.model.3)

fit_cis_95 <- confint(atl.model.3, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(atl.model.3, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.atl.3 <- bind_cols(coef.atl.3, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef
atl_boundary_white <- coef.atl.3[1,]


# austin
aus.model.3 = lm_robust(lfelonies_sd ~ aus_white_blv + aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                          hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                          lviolentcrime + phys_bound + com_density, 
                        data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)


aus.model.3_sense = lm(lfelonies_sd ~ aus_white_blv + aus_ses_blv + lpop + p_race_white + age_15_35_male + 
                         hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                         lviolentcrime + phys_bound + com_density, 
                       data = aus_fin, subset = pop > 0)


# create df of austin coefficients
coef.aus.3 <- tidy(aus.model.3)

fit_cis_95 <- confint(aus.model.3, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(aus.model.3, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.aus.3 <- bind_cols(coef.aus.3, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef
aus_boundary_white <- coef.aus.3[1,]


# boston
bos.model.3 = lm_robust(lfelonies_sd ~ bos_white_blv + bos_ses_blv + lpop + p_race_white + age_15_35_male + 
                          hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + 
                          phys_bound + com_density, 
                        data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)


bos.model.3_sense = lm(lfelonies_sd ~ bos_white_blv + bos_ses_blv + lpop + p_race_white + age_15_35_male + 
                         hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lcrime + 
                         phys_bound + com_density, 
                       data = bos_fin, subset = pop > 0)


# create df of bostin coefficients
coef.bos.3 <- tidy(bos.model.3)

fit_cis_95 <- confint(bos.model.3, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(bos.model.3, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.bos.3 <- bind_cols(coef.bos.3, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef
bos_boundary_white <- coef.bos.3[1,]

# chicago
chi.model.3 = lm_robust(lfelonies_sd ~ chi_white_blv + chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                          hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                          lviolentcrime + phys_bound + com_density, 
                        data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi.model.3_sense = lm(lfelonies_sd ~ chi_white_blv + chi_ses_blv + lpop + p_race_white + age_15_35_male + 
                         hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                         lviolentcrime + phys_bound + com_density, 
                       data = chi_fin, subset = pop > 0)

# create df of chicago coefficients
coef.chi.3 <- tidy(chi.model.3)

fit_cis_95 <- confint(chi.model.3, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(chi.model.3, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.chi.3 <- bind_cols(coef.chi.3, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef
chi_boundary_white <- coef.chi.3[1,]

# louisville - does not have misdemeanor vs felony distinction

# milwaukee
mil.model.3 = lm_robust(lfelonies_sd ~ mil_white_blv + mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                          hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                          lviolentcrime + phys_bound + com_density, 
                        data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

mil.model.3_sense = lm(lfelonies_sd ~ mil_white_blv + mil_ses_blv + lpop + p_race_white + age_15_35_male + 
                         hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                         lviolentcrime + phys_bound + com_density, 
                       data = mil_fin, subset = pop > 0)

# create df of milwaukee coefficients
coef.mil.3 <- tidy(mil.model.3)

fit_cis_95 <- confint(mil.model.3, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(mil.model.3, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.mil.3 <- bind_cols(coef.mil.3, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef
mil_boundary_white <- coef.mil.3[1,]


# seattle
sea.model.3 = lm_robust(lfelonies_sd ~ sea_white_blv + sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                          hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                          lviolentcrime + phys_bound + com_density, 
                        data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

sea.model.3_sense = lm(lfelonies_sd ~ sea_white_blv + sea_ses_blv + lpop + p_race_white + age_15_35_male + 
                         hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                         lviolentcrime + phys_bound + com_density, 
                       data = sea_fin, subset = pop > 0)


# create df of seattle coefficients
coef.sea.3 <- tidy(sea.model.3)

fit_cis_95 <- confint(sea.model.3, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(sea.model.3, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

coef.sea.3 <- bind_cols(coef.sea.3, 
                        fit_cis_95, 
                        fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")

# extract just white boundary coef
sea_boundary_white <- coef.sea.3[1,]


## combine 6 dfs (lousiville not in this analysis) w white boundary coefs
combined.coeffs <- rbind(atl_boundary_white, aus_boundary_white, bos_boundary_white, chi_boundary_white,
                         mil_boundary_white, sea_boundary_white)

combined.coeffs <- combined.coeffs %>% dplyr::select(-c(SE,statistic,p.value, conf.low, conf.high))

## make xaxis labels
Plot1labels <- c("Atlanta", "Austin", "Boston","Chicago","Milwaukee","Seattle")

## sensitivity analysis

# robustness value exercise 

sense_out1 = 
  sensemakr(atl.model.3_sense,
            treatment = "atl_white_blv",
            benchmark_covariates = "lviolentcrime",
            kd = 1)


sense_out2 = 
  sensemakr(aus.model.3_sense,
            treatment = "aus_white_blv",
            benchmark_covariates = "lviolentcrime",
            kd = 3)

sense_out3 = 
  sensemakr(bos.model.3_sense,
            treatment = "bos_white_blv",
            benchmark_covariates = "lcrime",
            kd = 1)

sense_out4 = 
  sensemakr(chi.model.3_sense,
            treatment = "chi_white_blv",
            benchmark_covariates = "lviolentcrime",
            kd = 1)

sense_out5 = 
  sensemakr(mil.model.3_sense,
            treatment = "mil_white_blv",
            benchmark_covariates = "lviolentcrime",
            kd = 1)

sense_out6 = 
  sensemakr(sea.model.3_sense,
            treatment = "sea_white_blv",
            benchmark_covariates = "lviolentcrime",
            kd = 1)


rvdf = data.frame(
  
  rv_val = c(as.numeric(sense_out1$sensitivity_stats$rv_q),
             as.numeric(sense_out2$sensitivity_stats$rv_q),
             as.numeric(sense_out3$sensitivity_stats$rv_q),
             as.numeric(sense_out4$sensitivity_stats$rv_q),
             as.numeric(sense_out5$sensitivity_stats$rv_q),
             as.numeric(sense_out6$sensitivity_stats$rv_q)),
  
  variable = combined.coeffs$Variable,
  
  bound = c("2x Violent \nCrime", "3x Violent \nCrime", "2x Crime", "415x Violent \nCrime", 
            "1793x Violent \nCrime","4x Violent \nCrime")
  
)


combined.coeffs$se = (combined.coeffs$conf.high_95 - combined.coeffs$Coefficient) / 1.96

out_meta = metagen(
  TE = combined.coeffs$Coefficient,
  seTE = combined.coeffs$se
)

# Atlanta: 2x violent crime
# Austin: 3x violent crime
# Boston: 2x crime
# Chicago: 415x violent crime
# Louisville: NA
# Milwaukee: 1793x violent crime
# Seattle: 4x violent crime


df_meta = data.frame(
  Variable = "Meta-analysis",
  Coefficient = out_meta$TE.random,
  df = out_meta$df.Q,
  outcome = "larrest_sd",
  se = out_meta$seTE.random,
  conf.low_95 = out_meta$TE.random - 1.96 * out_meta$seTE.random,
  conf.high_95 = out_meta$TE.random + 1.96 * out_meta$seTE.random,
  conf.low_90 = out_meta$TE.random - 1.645 * out_meta$seTE.random,
  conf.high_90 = out_meta$TE.random + 1.645 * out_meta$seTE.random
)

combined.coeffs = bind_rows(combined.coeffs, df_meta)
combined.coeffs$Variable = 
  factor(combined.coeffs$Variable, levels = c("Meta-analysis", "atl_white_blv","aus_white_blv","bos_white_blv",
                                              "chi_white_blv","mil_white_blv","sea_white_blv"))

fel_white_blv_plot_all_cities <- 
  ggplot(combined.coeffs, 
         aes(x = Variable, y = Coefficient)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 1) +
  geom_point(aes(x = Variable, 
                 y = Coefficient)) + 
  geom_errorbar(aes(x = Variable,
                    ymin = conf.low_90,
                    ymax = conf.high_90),
                width = 0,
                size = .8) +
  geom_errorbar(aes(x = Variable,
                    ymin = conf.low_95,
                    ymax = conf.high_95),
                width = 0,
                size = .4) +
  annotate("text",
           label =  paste0("RV: ", round(rvdf$rv_val[1:6], 2),
                           "\n", rvdf$bound[1:6]),
           y = combined.coeffs$Coefficient[1:6],
           x = c(2.6, 3.6, 4.5, 5.6, 6.6, 7.6),
           size = 3.75,
           family = "serif") +
  geom_hline(yintercept = combined.coeffs$Coefficient[7], 
             colour = gray(1/2), lty = 2,
             alpha = .4,
             size = .2) + ylim(-0.3, 0.8) + 
  ggtitle("Influence of White Racial Boundaries on Felony Arrests") + 
  xlab("Cities") + ylab("Racial Boundary Coefficients") + 
  scale_color_grey(start = .6, end = 0) + 
  theme_tufte(base_size = 14) + 
  scale_x_discrete(labels = c("Meta-Analysis",Plot1labels))


# make meta-analytic coefficient variable red
fel_white_blv_plot_all_cities <- fel_white_blv_plot_all_cities + 
  geom_point(data=subset(combined.coeffs, Variable == "Meta-analysis"),
             colour="red") +
  geom_errorbar(data = subset(combined.coeffs, Variable == "Meta-analysis"),
                mapping = aes(x = Variable,
                              ymin = conf.low_90,
                              ymax = conf.high_90),
                color="red",
                width = 0,
                size = 0.8) +
  geom_errorbar(data = subset(combined.coeffs, Variable == "Meta-analysis"),
                mapping = aes(x = Variable,
                              ymin = conf.low_95,
                              ymax = conf.high_95),
                color="red",
                width = 0,
                size = .4)

# expand size of plot so full labels can be seen
fel_white_blv_plot_all_cities <- fel_white_blv_plot_all_cities + 
  expand_limits(x= c(length(levels(combined.coeffs$Variable)) + 1))


ggsave(plot = fel_white_blv_plot_all_cities, filename = "fel_white_blv_plot_sens.png", 
       width = 8, height = 4, bg = 'white')


# put felony & misdemeanor plots side by side
library(gridExtra)
combined_misd_fel_plots = arrangeGrob(misd_white_blv_plot_all_cities, fel_white_blv_plot_all_cities, ncol = 2)


ggsave(plot = combined_misd_fel_plots, filename = "combined_misd_fel_plots.png", 
       width = 14, height = 6, bg = 'white')


